% Define approximate molar mass of components
molmass_casamino = 135.5; % g/L

% Total carbon molar mass
molmass_NZCasePlus = 125/molmass_casamino;
clc;
close all;
%% Establish concentrations to evaluate
ivec = [1:6]'; % Only use first 4 concentrations


%% 37C NZCasePlus
tempmat = monod_mat_NZCasePlus_25C;
smat = size(tempmat);

conc_25C = tempmat(ivec,1)/100*molmass_NZCasePlus*1000; % Convert to mM
grmax_25C = nanmean(tempmat(ivec,2:4),2);
grmax_err_25C = nanstd(tempmat(ivec,2:4),0,2);

%% 30C NZCasePlus
tempmat = monod_mat_NZCasePlus_30C;
smat = size(tempmat);

conc_30C = tempmat(ivec,1)/100*molmass_NZCasePlus*1000; % Convert to mM
grmax_30C = nanmean(tempmat(ivec,2:4),2);
grmax_err_30C = nanstd(tempmat(ivec,2:4),0,2);

%% 37C NZCasePlus
tempmat = monod_mat_NZCasePlus_37C;
smat = size(tempmat);

conc_37C = tempmat(ivec,1)/100*molmass_NZCasePlus*1000; % Convert to mM
grmax_37C = nanmean(tempmat(ivec,2:4),2);
grmax_err_37C = nanstd(tempmat(ivec,2:4),0,2);

%% Now do Arrhenius analysis

conc_tot = conc_37C; % All should be the same concentration
T = [25 30 37]';

grmax_tot_NZCasePlus = [grmax_25C grmax_30C grmax_37C];
grmax_err_tot_NZCasePlus = [grmax_err_25C grmax_err_30C grmax_err_37C];

figure;

Ea_NZCasePlus = [];
Ea_err_NZCasePlus = [];

for k=1:length(ivec)
    g = grmax_tot_NZCasePlus(k,:)';
    ge = grmax_err_tot_NZCasePlus(k,:)';
    [f1,f2] = fit(1./(T+273.15), log(g),'poly1');
    errorbar(1./(T+273.15), log(g), ge./g, 'ko', 'linewidth', 1);
    hold on;
    plot(1./(T+273.15),f1(1./(T+273.15)), 'k', 'linewidth', 1);
    
    disp('Activation Energy (KT):');
    disp(-f1.p1/298);
    
    disp('With error:');
    
    alpha = 0.95;
    ci = confint(f1, alpha);
    fiterr = abs(f1.p1 - ci(1,1))/298; % In units kbT
    disp(fiterr);
    
    Ea_NZCasePlus = [Ea_NZCasePlus; -f1.p1/298];
    Ea_err_NZCasePlus = [Ea_err_NZCasePlus; fiterr];
    
end

set(gca, 'fontsize', 20);
set(gcf, 'Position', [0 0 400 300]);
xlabel('1/T (1/K)');
ylabel('Log (growth rate)');
xlim([1/312 1/296]);

